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ABSTRACT 

Time-dependent, axisymmetric hydrodynamic simulations have been used to study accretion disks 
consisting of counterrotating components with an intervening shear layer (s). Configurations of this type 
can arise from the accretion of newly supplied counterrotating matter onto an existing corotating disk. 
The grid-dependent numerical viscosity of our hydro code is used to simulate the influence of a turbulent 
viscosity of the disk. Firstly, we consider the case where the gas well above the disk midplane (z > 0) 
rotates with angular rate -1-17 (r) and that well below (z < 0) has the same properties but rotates with 
rate —VL(r). We find that there is angular momentum annihilation in a narrow equatorial boundary layer 
in which matter accretes supersonically with a velocity which approaches the free-fall velocity. This is 
in accord with the analytic model of Lovelace and Chou (1996). The average accretion speed of the 
disk can be enormously larger than that for a conventional a— disk rotating in one direction. Under 
some conditions the interface between the corotating and counterrotating components shows significant 
warping. Secondly, we consider the case of a corotating accretion disk for r < rt and a counterrotating 
disk for r > vt- In this case we observed, that matter from the annihilation layer lost its stability and 
propagated inward pushing matter of inner regions of the disk to accrete. Thirdly, we investigated the 
case where counterrotating matter inflowing from large radial distances encounters an existing corotating 
disk. Friction between the inflowing matter and the existing disk is found to lead to fast boundary layer 
accretion along the disk surfaces and to enhanced accretion in the main disk. 

These models are pertinent to the formation of counterrotating disks in galaxies and possibly in Active 
Calactic Nuclei and in X-ray pulsars in binary systems. For galaxies the high accretion speed allows 
counterrotating gas to be transported into the central regions of a galaxy in a time much less than the 
Hubble time. 

Subject headings: accretion, accretion disks — galaxies: formation — galaxies: evolution — galaxies: 
nuclei — galaxies: spiral — galaxies 



1. INTRODUCTION 

The widely considered models of accretion disks have gas 
rotating in one direction with a turbulent viscosity acting 
to transport angular momentum outward (Shakura 1973; 
Shakura & Sunyaev 1973). However, recent observations 
indicate that there may be more complicated disk struc- 
tures in some cases on a galactic scale, in galactic nuclei, 
and in disks around compact stellar mass objects. Stud- 
ies of normal galaxies have revealed counterrotating gas 
and/or stars in many galaxies of all morphological types 
- ellipticals, spirals, and irregulars (Rubin 1994a, 1994b; 
Galletta 1995). In elliptical galaxies, the counterrotating 
component is usually in the nuclear core and may result 
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from merging of galaxies with opposite angular momen- 
tum. 

In a number of spirals and SO galaxies, counterrotating 
disks of stars and/or gas have been found to co-exist with 
the primary disk out to large distances (10 — 20 kpc), with 
the first example, NGC 4550, discovered by Rubin, Gra- 
ham, & Kenney (1992). Of interest here is the "Evil Eye" 
galaxy NGC 4826 in which the direction of rotation of the 
gas reverses going from the inner (180 km/s) to the outer 
disk (—200 km/s) with an inward radial accretion speed 
of more than 100 km/s in the transition zone, while the 
stars at all radii rotate in the same direction as the gas 
in the inner disk which has a radius ~ 1200 pc (Braun et 
al. 1994; Rubin 1994a, 1994b). The large scale counter- 
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Fig. 1. — Sketches of three main cases considered in this paper. The signs ® indicates corotating matter and 
counterrotating matter. Details of the models are described in the text. 



rotating disks probably do not result from mergers of flat 
galaxies with opposite spins because of the large vertical 
thickening observed in simulation studies of such mergers 
(Hernquist & Barnes 1991; Barnes 1992). Thakar & Ryden 
(1996) discuss different possibilities, (a) that the counter- 
rotating matter comes from the merger of an oppositely 
rotating gas rich dwarf galaxy with an existing spiral, and 
(b) that the accretion of gas occurs over the lifetime of the 
galaxy with the more recently accreted gas counterrotat- 
ing. Subsequent star formation in the counterrotating gas 
may then lead to counterrotating stars. Counterrotating 
gas may have a strong interaction with the gas resulting 
from mass loss by corotating stars (Braun ct al. 1994). 
Further, the two-stream instability between counterrotat- 
ing gas and corotating stars can be important for gas ac- 
cretion (Lovelace, Jore, & Haynes 1996). 

Accretion of counterrotating matter onto an existing 
corotating disk around a massive black hole may occur in 
an Active Galactic Nucleus (AGN). The counterrotating 
matter may come from recently captured material from a 
galaxy merger, from a tidally disrupted molecular cloud, 
or from a star. Tidal disruption of a star by the central 
black hole may release essential part of the star in the form 
of gas (Rees 1990; Shlosman, Begelman, & Frank 1990; 
Khokhlov & Melia 1996). There is a significant probabil- 
ity that this gas has angular momentum opposite to that 
of the main disk. This can lead to formation of small disks 
with opposite angular momentum in the center of a galaxy 
(Melia 1994). This was studied using 3D numerical sim- 
ulations by Ruffert (1994), Ruffert & Melia (1994). The 
retrograde disks may also slow down the rotation of the 
central black hole (Moderski & Sikora 1996). The greatly 
enhanced accretion rate for counterrotating matter in the 
presence of a corotating disk may give outbursts of AGN. 

Another situation where counterrotating matter may 
encounter an existing corotating disk is in low mass X-ray 
binary sources where the accreting, magnetized rotating 
neutron stars are observed to jump between states where 
they spin-up and those where they spin-down. Nelson et 
al. (1997a, b) and Chakrabarty et al. (1997) have pro- 
posed that the change from spin-up to spin-down results 
from a reversal of the angular momentum of the wind sup- 
plied accreting matter. Their suggestion is based on the 



recent BATSE observations which show a rapid change 
from spin-up to spin-down evolution of some X-ray pulsars. 
This is difficult to explain with conventional theory (see, 
however, Lovelace, Romanova, & Bisnovatyi-Kogan 1998). 
Reversal of the angular momentum of accreted matter 
in wind-fed binaries was observed in the two-dimensional 
simulations by Blondin et al. (1990), Livio et al. (1991), 
and Bisikalo ct al. (1996). 

An important open problem is to understand the 
hydrodynamic interaction of counterrotating gaseous disks. 
Earlier two-dimensional simulations used either collision- 
less particles to represent a system of stars (Davies & 
Hunter 1997) or "sticky" particles (Thakar & Ryden 1996, 
and references therein) to represent a system of clouds. 
The sticky particle interaction gives only a rough repre- 
sentation of a viscous fluid. When both components of 
a counterrotating system are gaseous, it is important to 
understand their interaction using the full equations of 
hydrodynamics including the viscous transport of angular 
momentum. Hydrodynamic interaction of a galactic disk 
with infalling matter of a non-rotating corona was inves- 
tigated by Falcke & Melia (1997) using a one dimensional 
approach. They showed that this interaction may lead to 
significant evolution of the disk. 

Lovelace & Chou (1996) (hereafter LC) developed an 
analytic theory of the viscous hydrodynamic interaction 
of counterrotating gaseous disks basing on the a— viscosity 
model of Shakura (1973) and Shakura & Sunyaev (1973). 

Here, we present the results of two-dimensional, time- 
dependent hydrodynamic simulations of counterrotating 
accretion disks. We investigate three main models which 
are sketched in Figure 1. 

Model I is shown in the top sketch of Figure 1: In this 
case the accretion disk consists of counterrotating gaseous 
components with gas above the disk midplane rotating 
with angular rate +Q{r), and that below rotating at rate 
— fi(r). The interface between the components at 2; ~ 
constitutes a supersonic shear layer. A configuration of 
this type may arise from the accretion of newly supplied 
counterrotating gas onto an existing corotating disk. For 
this case the analytical model of Lovelace & Chou (1996) 
proposed that matter approaching the equatorial plane 
from above and below has angular momenta of opposite 
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signs with the rcsuh that there is angular momentum anni- 
hilation at z ~ 0. This matter loses its centrifugal support 
and accretes at essentially free-fall speed. 

Model II is shown in the middle sketch in Figure 1 : In 
this model we consider that where there is a corotating 
(+ri) accretion disk (the "old" disk) for r < rt and a 
counterrotating (— fi) disk (the "new" disk) in the region 
r > rt- This case is clearly pertinent to the "Evil Eye" 
galaxy NGC 4826. 

Model III is shown in the lower sketch of Figure 1: 
This model corresponds to inflow of counterrotating mat- 
ter which interacts with an existing corotating disk. In 
this case annihilation of angular momentum occurs on the 
surface of the existing disk and this can lead to fast accre- 
tion. 

In §2 we discuss the pertinent fluid dynamical equations 
for axisymmetric accretion, the initial and boundary con- 
ditions, and the numerical method used. In §3 - §5 we 
discuss results obtained for Models I-III. We present as- 
trophysical examples in §6. In §7 we give conclusions of 
this work. 



Here, p{R,0) is the density; v{R,9) the flow veloc- 
ity, with components {vr,v^,V0)\ p{R,6) the pressure; 
$g the gravitational potential; £ the full speciflc energy 
£ = e + {v% + v^ + V0)/2; e the specific internal energy; and 
w the full specific enthalpy, w = e+p/p+{v'jf^ + v^ + Vg)/2. 
An equation of state is necessary, and we assume the ideal 
gas equation, 

p= (7- l)ep, 

where 7 is the usual ratio of heat capacities. The gravita- 
tional potential is considered to be the Newtonian poten- 
tial of the central object. 

GM 



That is, we assume that the self gravity of the disk is neg- 
ligible. 



2. MATHEMATICAL MODEL 

2.1. Basic Equations 

The axisymmetric flows are described by the inviscid hy- 
drodynamic equations (Euler equations) in spherical, iner- 
tial coordinates {R, (j), 9) (with ^ = the equatorial plane): 
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2.2. Alternative Form of Energy Equation 

In our early disk simulations we found that the energy 
equation taken in divergent form (le) can lead to inac- 
curate results. The reason for this is the large value of 
the azimuthal velocity and the corresponding kinetic 
energy in comparison with other velocities/energies. 
The relatively small quantity e is determined from (le) 
inaccurately in this case. 

To correct this problem we now have excluded the (/>- 
velocity part of the energy from the energy equation by 
writing 
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and subtracting this equation from equation (le) to obtain 
another divergent form of energy equation. 
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where £' and u/ are the 'reduced' (poloidal) full specific 
energy and the 'reduced' (poloidal) full specific enthalpy, 



£' 



€+{vji + vj)/2 and iv' = e + p/p+ {v% + v^)/2, re- 



spectively. In all of our calculations equation (2) was used 
(le) instead (le). 
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Fig. 2. — Log-scaled iso-density lines from pmin = 10~^ to Pmax = 1 for a Keplerian disk described in §2.4.1 with 
CgQ = 0.01 and 7 = 1.01. The dashed lines show the standard disk scale-height, H = tCs/Vk- The circles represent the 
region + z'^ < Rq where matter is considered to be absorbed. 
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Fig. 3. — Log-scaled iso-density lines from Pmin — 

10-6 to Pmax = 1 for a sub-Keplerian disk as described in §2.4.2 for 

A = 1/3, n = 11.5i?o, and n = 3/2 (or 7 = 5/3). 
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Fig. 4.- Example of nature of grid used in simulations. The grid is in spherical coordinates with R the radial distance 
and the co-latitude measured from the equatorial plane. The grid is non-homogeneous with increments increasing 
exponentially with both R and |^|. The grid shown has Nji x Ng = 50 x 50. 



2.3. Characteristic Scales 

The physical quantities have evident characteristic 
scales: Distances are measured in terms of the inner ra- 
dius of the disk Rq. Velocities arc measured in terms 
of the Keplerian velocity at the inner edge of the disk, 
Vko = \/ GM/ Ro . Time is measured in units of the period 
of revolution of the inner edge of the disk to = 2'kRq/Vkq- 
The density is given in relative units. 

2.4. Initial Equilibria 

We have adopted spherical coordinates for our simu- 
lations in order to obtain a good match of the grid to 
the accretion disk which in general has a thickness in the 
2— direction {2H) increasing with R. However, for conve- 
nience our figures and initial disk equilibria described in 
this subsection are given in (r, z) coordinates. 



A disk equilibrium with v. 
described by the equations 



2.4.1. InitAal Conditions 1 

0, Vz = 0, but 



1 dp 
p dr 

1 dp 
p dz 



r 



dr 



dz 



(3a) 
(3b) 



where the position and velocity are now dimensionless as 
mentioned. 

With the assumption of a polytropic gas p = Kp^^^/'^, 
where n is the polytropic index, n = 1/(7—1), the common 
solution of equations (3) can be written as 



w • 



Wo + F - ^„ = Wo + F + 



;2 ' 



(4) 



where w is the specific enthalpy, w = K{n+l)p^^"' = nc^, 
with wq a constant, Cs — \/lvl P the sound speed, and 
F(r) an arbitrary function depending only on r. The an- 
gular velocity is given in this case as 



v'i = r 



dF 
dr 



The disk boundary is determined by the condition w = 
or p = 0. If we assume the Keplerian law for angular 
velocity, t;? = 1/r, then 



1 



= $g(r,^ = 0). 



and 



P = Po 



1- 



(5) 
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where po and CsQ are the values at the point {r = 1, z = 0). 
The disk boundary in this case is 



1 



SO' 



so the disk thickness diverges at some very large radius 
outside the region of interest. A sample density distribu- 
tion is shown in Figure 2. Hereafter, we refer to this case as 
Initial Condition 1. This was used for the runs described 
in §3 and §4. 

2.4.2. Initial Conditions 2 

A second equilibrium disk model was used in our sim- 
ulations. For this we took the disk boundary (p = 0) to 
have the form 



1 



where ri and A are parameter determining disk radial size 
and thickness. The rotation (slightly sub-Keplerian), en- 
thalpy, and density for this case are: 



1- V 



3/2 



W 



P = Pqw 



= 



A2 



1 + A2 ■ 



(6a) 



(6b) 



(6c) 



The density distribution for this case is shown in Figure 
3. Hereafter, we refer to this case as Initial Conditions 2. 

For building equilibria we have used both a rotation 
law based approach and one based on the shape of the 
disk smface. For the first case, the equilibria are found 
for known rotation law according to the chain v^{r) i— * 
F{r) I— > w{r,z) i— > z{r)\p^Q. For the second case, a known 
disk surface is used according to the chain z{r)\p^Q i—^ 

F{r) = ^{r,z{r)) — wo ^^l'')- Abakumov et al. (1996) 
use an analogous approach. 



2.6. Numerical Method and Viscosity Treatment 

An explicit Lax-Fricdrichs finitc-diffcrcncc scheme with 
flux correction in Chakravarthy-Osher (Chakravarthy & 
Osher 1985) form was used. This scheme is 1st order of 
approximation in time and 3rd order of approximation in 
space and is oscillation-free near discontinuities. The com- 
bined Lax-Friedrichs-Osher scheme was suggested in work 
by Vyaznikov, Tishkin, & Favorsky (1989). 

Figure 4 shows the non- homogeneous {R, 0) grid used 
in the simulations. The radial grid increment 5R increases 
outwards exponentially, SR = SR^a^^^", where a = 1.05. 
The 6 grid 59 increases in a similar way going away from 
the equatorial plane. For a disk which has thickness h 
increasing with r, this type of grid gives a more uniform 
distribution of cells within the disk including its central 
regions. In the simulations we used grids of A^j? x Ng of 
100 X 100 or 200 x 200 in some cases. 

The simulation method used naturally has numerical 
viscosity. With a fixed grid, we cannot vary the numer- 
ical viscosity, but we can investigate its influence. The 
coefficient of numerical viscosity can be represented as 
i^num ~ CgAh, where Cg is the sound speed. Ah is the size 
of the grid. This expression for v is analogous to that pro- 
posed by Shakura (1973) and Shakura & Sunyaev (1973) 
for the turbulent viscosity of accretion disks, v = aCgH , 
where H is the thickness of the disk, and a^lO^^ — 10^^ 
is a dimensionless parameter. We calculated the steady 
accretion for a disk rotating in one direction and observed 
the accretion due to the numerical viscosity. Using the 
formula for the accretion speed of a viscous accretion disk, 
Vr ~ y/r, we estimated that a < 1 for r > 5i?o for a 
200 X 200 grid and a < 1 for r > 7J?o for a 100 x 100 
grid, and a decreases to a ^ 0.03 — 0.1 at the outer radii. 
A challenge for future work is the inclusion of a true a 
viscosity. 

We investigated three different models of counterrotat- 
ing accretion flows which correspond to different possible 
astrophysical situations. The first model corresponds to 
disks counterrotating in the z— direction (Model I, §3). 
The second corresponds to disks counterrotating in the 
r— direction (Model II, §4). The third corresponds to coun- 
terrotating gas infalling onto a pre-existing corotating disk 
(Model HI, §5). 



2.5. Boundary Conditions 
Simulations were performed in the region Rq < R < 

Rmax, -7r/2 <e< 7^/2. 

The outer boundary conditions R = Rmax are treated 
following the method used by Sawada based on solving 
of the Riemann problem (see, e.g., Sawada, Matsuda, & 
Hachisu 1986; Sawada & Matsuda 1992). These bound- 
ary conditions were used for Models I and II. However, 
for Model HI different boundary conditions are applied in 
that matter is considered to inflow supersonically through 
part of the outer boundary. 

We assume that all matter coming into the sphere 
R ^ -Ro is absorbed. On the z— axis, = ±7r/2, we of 
course have vo = and = 0. 

No symmetry is assumed with respect to the equatorial 
plane. 



3. DISK COUNTERROTATING IN THE Z-DIRECTION: 

MODEL I 

Here, we discuss simulation results for Model I where 
initially the upper half of the disk (z > 0) is corotating 
{+Qk) while the lower half {z < 0) is counterrotating 
{—V,k)- This case was treated analytically by LC, and 
therefore we have a test both of the LC theory and of our 
simulation code. 

First in §3.1, we check the initial equilibrium disk for 
the conventional case of unidirectional disk rotation. Sec- 
ondly in §3.2, we investigate the counterrotating disk for 
7 = 1.01, wliic;li is close to isothcirmal. This value of 7 close 
to unity mimics conditions of strong radiative cooling. In 
this subsection we compare the simulations results with 
the predictions of LC where radiative cooling is included. 
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Fig. 5. — Model I. Iso-density lines for a conventional disk rotating in one direction for Initial Conditions 1 after a time 
t = 40to, where to is the period of rotation of the inner edge of the disk at r = i?o- The disk is close to its initial 
configuration, excluding the inner part, which has changed as a result of numerical viscosity. Simulations were performed 
with a grid 100 x 100. 
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Fig. 6. — Model I. Iso-density contours (solid lines) and mass flux vectors (pv) in counterrotating accretion disk at time 
t = 14to, where to is the rotation period at the inner radius of the disk. The poloidal flow near the midplane of the disk 
is supersonic, excluding the small regions separated by the bold lines. Simulations were performed with a grid 100 x 100. 
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r/Ro 

Fig. 7. — Model I. Variation of radial velocity Vr with r in the midplane of the disk at times t = 6to, 14fo, and 32^0, 
where to is the rotation period at the inner radius of the disk. 




r/Ro 

Fig. 8. — Model I. Variation of density with r in the midplane of the disk at the same time intervals as Figure 7. 



In these subsection we used a 100 x 100 grid. Thirdly in 
§3.3, we discuss the cases 7 = 1.1 and 5/3 which corre- 
sponds to non-isothermal conditions where numerical vis- 
cous heating is significant. In this subsection we used a 
200 X 200 grid. 

3.1. Conventional Accretion Disk 

As a calibration of our simulation code, wc investigated 
a conventional ac;cretion disk which rotates in only one di- 
rection. The disk evolves as a result of the relatively small 
numerical viscosity of the code. We took Initial Conditions 
1 (see §2.4.1) with 7 = 1.01 and performed long-term nu- 



merical simulations of this disk. We observed that the disk 
shape did not change during many periods of rotation of 
the inner radius of the disk. Figure 5 shows the iso-density 
contours distribution at t = 40io- The evolved disk is close 
to that at i = (sec Figure 2), excluding the inner part, 
where the thickness of the disk is very small (compara- 
ble with the grid size), and the numerical viscosity has a 
noticeable influence. 

3.2. 7 = 1.01 

This case is close to isothermal and is closest to the LC 
analytical model where both viscous heating and radiative 
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Fig. 9. — Model I. Radial velocity (hollow circles) and azimuthal velocity (filled circles) velocities across the disk 
as a function oi ( = z/h at radial distance r = TRq. Here, h{r) is the characteristic half-thickness of the region of the 
strongest inflow, and = GM/r. The solid lines show the theoretical dependencies from LC. Note that to a good 
approximation Vr is an even function of C and is an odd function. 




Fig. 10. — Model I. Mass accretion rate as a function of radial distance r at different moments of time. The upper curves 
are for a counterrotating disk at times t — Mq, 7to, and 14fo. The lower curve (C) is for a conventional accretion disk 
rotating in one direction at time t = 40to, where the accretion results from the numerical viscosity of the code. 




Fig. 11. — Model I. Contours of density and velocity vectors for 7 = 1.1 at times t = 4.8to (bottom panel) and t = 6.2to 
(top panel). The bold line separates the corotating and counterrotating gas. The simulations were performed for the part 
of the disk 3 < r/i?o < 7 with a 200 x 200 grid. 



cooling are included. 

We took Initial Conditions 1 for the initial conditions for 
the density and velocity of the disk (see Figure 2), but put 
the angular velocity in the lower-half of the disk (z < 0) 
opposite to the angular velocity of the upper-half {z > 0). 
Figure 6 shows that matter starts to rapidly flow radially 
inward in a narrow region near the midplane of the disk. 
The poloidal flow near the midplane is supersonic, and the 



radial inflow velocity increases as r decreases. 

Figure 7 shows the radial variation of the radial velocity 
in the midplane of the disk at different times. The ve- 
locity distribution changes only gradually with time. The 
velocity close to the gravitating center, r/Ro ^ 1 — 3 is 
much larger than velocity in the outer regions of the disk, 
r/i?o ~ 16 — 18. Note, that for r ~ i?o, the radial inflow 
velocity is close to the free-fall velocity V// = \/2GM/R. 
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Fig. 12. — Model I. Profiles of the density (1), pressure (2), and temperature (3) through the disk at r = 5Ro at the time 
t = 6.2to corresponding to the top panel of Figure 11. 



An inflow velocity ^ GM/R is predicted by LC. Veloc- 
ity is strongly supersonic with Mach number M > 100. 
Figure 7 shows that velocity distribution along the disk is 
almost the same during the long time. Figure 8 shows the 
dependence of density on radial distance. 

The density of the disk decreases with time, specifically 
in the inner regions of the disk. Both, the velocity and 
density show small-scale wave-like variations as a function 
of radius. 

Figure 9 shows the observed vertical (z) variation of the 
radial and azimuthal velocities in the disk at r = TRq and 
the theoretical dependencies calculated by LC. The dimen- 
sionless vertical length was introduced by LC as ( = z/h, 
where h = h{r) is the characteristic half-thickness of the 
boundary layer in which most of the radial inflow of mat- 
ter occurs. Consistent with LC, we find that {h/H)"^ <C 1 
over most of the disk surface, where H is the full disk 
half-thickness. At r = 7Ro, we find h/H « 0.25. One can 
sec that radial inflow velocity has a maximum in the mid- 
plane of the disk and decreases to small values at large 
Note that we observe weak outflows for 4 < ^ < 8 where 
Vr > 0. Similar outflows were predicted by LC (the solid 
curves). The azimuthal velocity grows with increasing 

( and reaches the Keplerian velocity for C ~ 2.4. There 
is also a region, where the azimuthal velocity is slightly 
larger than Keplerian velocity, 3 < C < 4.5. A similar 
region was predicted by LC (sec solid curve in Figure 9). 

The conditions considered by LC were different from the 
initial conditions considered for our simulations. Thus, it 
appears that the main features of vertically stratified coun- 
terrotating flows do not depend strongly on the initial con- 
ditions. At longer times in our simulation runs, the region 
of fast equatorial flow becomes somewhat thicker. Over- 
all the simulation results are close to the predictions of 
the analytic model of LC. Although the analytic model in- 
cludes viscous heating and radiative cooling, these almost 
compensate each other with the result the model is close 
to the simulations which has almost isothermal conditions 
and no radiative cooling. 

Figure 10 shows the radial dependence of the mass ac- 
cretion rate M(r) at different times, and that for a disk 
rotating in one direction at time t = 40to (§3.1). The mass 
accretion rate for the counterrotating disk can be written 
as 

M ^ 2Tirp{r.z = Q)VK{r)[i.2hh{r)] , (7) 

where the numerical coefficient (3.25) is obtained from 
Figure 9. In that the surface mass density of the disk 
is S = 27rp(r, z = 0)2H, the average accretion speed is 



wcfi ~ l-Q2(h/H)VK- For comparison the mass accretion 
rate for a conventional a— disk rotating in one direction 
(Shakura 1973; Shakura & Sunyaev 1973) is 

Mss « 2nrp{r, z = 0)^^(0 (^^^ 2H(r) . (8) 

In this case the average accretion speed is uss = 
a{cs/VK)'^VK ■ For a conventional accretion disk, the 
quantity q;(cs/Vr-)^ is commonly thought to be much 
smaller than unity. For example, for Cs/Vk = 0.01, 
h/H = 0.1, and a = 0.1, M/Mss « 1-6 x 10^. 

3.3. 7 = 1.1 and 7 — 5/3 

We also did simulations of vertically stratified counter- 
rotating disks for 7 = 1.1 and 7 = 5/3 which correspond to 
non-isothermal conditions were the affect of viscous heat- 
ing is more important. We find in general that after a given 
time and for a given numerical grid, the temperatures in 
the disk increase as 7 — 1 increases. Thus, the disk thick- 
ness H/r is larger for larger 7 — 1. On the other hand for 
fixed 7, the temperature in the disk decreases gradually 
as the resolution x Ng increases. Here, we performed 
simulations with the grid 200 x 200. To increase the res- 
olution even more (to decrease viscosity) we performed 
simulations of only part of the disk: 3i?o <r < TRq. 

At early times (i/io ~ 1 — 3), the flow in the equa- 
torial plane of the disk was similar to that for 7 = 1.01. 
However, for longer times the region of fast poloidal inflow 
begins to warp. Figure 11 shows contours of density and 
velocity vectors of the disk at two times. The amplitude 
of the warp of the interface Az between corotating and 
counterrotating flows grows. After a fixed time the warp 
amplitude is a strongly increasing function of 7 — 1. For 
7 = 1.01 the fractional warp amplitude A^/r is very small 
for all times. However, the warp amplitude for 7 = 1.1 is 
Az/r « 0.034, and for 7 = 5/3 it is Az/r « 0.13 both at 
t/tQ = 7 with the same resolution as used in Figure 11. 
The warp amplitude Az is of the order of the disk half- 
thickness H. For 7 = 1.1 the dominant radial wavelength 
is Xr ~ 1.5i?o while for 7 = 5/3 there is a wide spec- 
trum of wavelengths from ^ Rq to ^ Rq. The warping is 
suppressed by low grid resolution. 

The warping may represent a Kelvin-Helmholtz type of 
instability driven by the shear in the poloidal fiow. The 
tidal gravitational force of the central object, F^, acts to 
keep the warp amplitude from growing without bound. 
Note that Kelvin-Helmholtz instabilities driven by the 
shear in the azimuthal flow are suppressed in our sim- 
ulations which are axisymmetric. Alternatively, the disk 
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warping may be duo to misymmctrical (top/bottom) heat- 
ing of the disk which gives an excess pressure in say the 
bottom half of the disk which displaces the disk upward. 
Plots of the pressure, density, and temperature through 
the disk as shown in Figure 12 support this simple expla- 
nation. 

4. DISK COUNTERROTATING IN THE R-DIRECTION: 
MODEL II 

Here, we treat Model II which is shown schematically 
in the middle panel of Figure 1. In this model the inner 
part of the Keplerian disk (r < r*) rotates in the positive 
direction whereas the outer part (r > rt) rotates in the 
opposite direction. The transition between positive and 
negative azimuthal velocity was centered at r/Ro = 14 
and extended over the narrow region 13.5 < r/Ro < 14.5. 
This transition region is small, involving only 5 grid points 
in the radial direction for a mesh size of 100 x 100. The 
boundary condition at the right side of the disk was taken 
the same as in Model I. Simulations were performed for 
7 = 1.1 and 7 = 5/3. Results are presented for 7 = 1.1. 

In our simulations, matter in the transition zone with 
low angular momentum flows rapidly inward compressing 
matter at smaller radii as shown in Figures 13 and 14. This 
inflow sets up a density wave which propagates inward in 
the disk. The matter behind the wave has lower density. 
The speed of propagation of the wave grows with time 
and reaches a maximum value equal to about one half of 
the free-fall velocity (see Figure 14). The azimuthal veloc- 
ity in the wave is sub-Keplerian with ^ (0.6 — 0.8)vk 
in inner part of the disk. Matter behind the wave (at 
larger r) has opposite velocity, with a magnitude some- 
what smaller than the local Keplerian velocity (sec Figure 
14). The wave propagates supersonically, with Mach num- 
ber M ~ 20 — 30. The Mach number grows as the wave 
moves inward (see Figure 14). 

The shock wave propagating initially at times t < lOto 
is not a standard shock wave, because it has significant 
angular momentum, and its formation and propagation is 
based on the deficit of angular momentum of the matter of 
the wave and an inbalancc of gravitational and centrifugal 
forces. The density and velocity jumps are much larger 
than in the case of standard shock wave (see Figure 14). 
This wave resembles a soliton-type wave. A similar type of 
wave was observed in accretion disks where a local deficit 
of angular momentum was caused by magnetically driven 
outflows (Lovelace, Romanova, & Newman 1994; Lovelace, 
Newman, & Romanova 1997). Note that the solid line in 
Figure 13 indicates the boundary between corotating and 
counterrotating components. The density wave reaches 
r = 8R0 during t = lOto (ten periods of rotation of the 
inner edge of the disk). For longer times, matter started 
to "reflect" from the inner regions and move outward. 

The effective collapse of the inner disk can be under- 
stood by considering the details of the model used. The 
density of the disk is constant along the equatorial plane 
so that the mass of the equidistant rings AR = const in- 
creases outwards, AM = 27rprAr ^ r, so that the mass 
of ring at r = 14 is 14 times larger than a ring at r = 1. 
At the same time the angular momentum of the rings in- 
creases outwards as AL AMvKr ^ r^/^, so that the 
angular momentum of inner ring is 52 times smaller than 



of the ring at r = 14. The layer between counterrotating 
disks has pretty large mass and zero angular momentum. 
It is not centrifugally supported and falls down to the cen- 
ter. It accumiilates the matter of the inner regions of the 
disk and their angular momentum. But absorbed angular 
momentum is not sufficient to rich the centrifugal barrier. 

The thickness of the boundary layer in which there is an- 
nihilation of angular momentum depends on the viscosity. 
For a larger viscosity, the layer is thicker and the propa- 
gating wave is stronger. The fast inward propagation of 
such a wave in a disk around a black hole would give an 
outburst in the luminosity as the wave reaches the black 
hole. This represents a possible mechanism of generating 
outbursts in Active Galactic Nuclei, for example. 

5. INFLOW OF COUNTERROTATING GAS ONTO A 
COROTATING DISK: MODEL III 

Here, we discuss Model HI which is shown schematically 
in the bottom panel of Figure 1. In this case counterro- 
tating gas of relatively low density inflows with signiflcant 
poloidal speed and encounters an existing Keplerian disk. 
For the simulations, we took Initial Conditions 2 (§2.4.2) 
corresponding to a finite radius, slightly sub-Keplerian, 
corotating disk (see Figure 3). The effective outer radius 
of this disk was taken to be at rout = 11.5i?o. At the 
outer, right-hand boundary of the computational region, 
matter is launched into the simulation region from the 
part of the boundary — 3 < z/Rq < 3 with radial veloc- 
ity Vr = 0.7Vk and azimuthal velocity = —Vk- When 
both: the "target" disk and incoming matter rotate in the 
same direction, then incoming matter moves inward only 
until it reaches the centrifugal barrier (see white line on 
Figure 15), then it stops, and moves in the ±z directions 
away from the disk surface. 

Figure 15 shows the observed behavior of density and 
velocity in this case. Note, that in this case we observe 
only very small accretion connected with finite viscosity 
of the disk. This case was considered for reference. 

When the inflowing matter rotates in the opposite di- 
rection to that of the "target" disk, then it also reaches 
a centrifugal barrica'. However, a non-negligible fraction / 
of the inflowing matter interacts with the disk causing an- 
gular momentum annihilation. For this case the fraction 
is / ~ 0.2. Figures 16 and 17 show the evolution. One 
can see from Figure 16 (top and bottom panels), that at 
times 2> < t/t{) < 7, the gas starts to flow along the sur- 
faces of the target disk. Later, at 7.5 < t/to < 15, it flows 
strongly along the surfaces (see top panel of Figure 17). 
The radial velocity of matter flowing along the surfaces is 
^ (0.3 — 0.5)Vff. Later, the layer became thicker owing 
to viscous dissipation and the flow became more chaotic 
(Figure 17, bottom panel). The temperatures increased 
few times in the surface layer. Here, we should remind, 
that viscous heating is included, but cooling is not. If ra- 
diative cooling balances viscous heating, then this surface 
layer may survive much longer than in the simulation. The 
surface flow is similar in some respects to that observed in 
Model I (§3). As in Model I, we also observed the increase 
of the radial inflow velocity with decreasing r. 

In this model, we also observed an interaction similar 
to that investigated in Model II. Namely, the interaction 
occurred not only along the surfaces of the disk but also 
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Fic. 13. Model II. Gray-scale plots of density and matter flux vectors showing the evolution of Model 11. Starting from 
the top, the panels show times t = 2tQ, AIq, QIq, 8to, and lOfo- The vectors are proportional to the value of the mass flux 
pv. The black solid line shows the line of zero azimuthal velocity. For the case shown, 7 = 1.1. 
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Fig. 14. — Model II. Radial distributions of density (top left panel), radial velocity (in units of local free-fall velocity) 
(top right panel), azimuthal velocity (in units of local Keplerian velocity) (bottom left panel), and Mach number (right 
bottom panel) at different times, t = 2to, 4to, 6^0, 8io, and lO^o- The time sequence corresponds to that shown at Figure 
13. 
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Fig. 15. — Model III. Results of simulations of Model III where radially inflowing matter encounters an existing corotating 
disk. The inflowing matter rotates in the same direction as the disk and there is no accretion. The inflowing matter is 
injected at 2 ~ ±3-Ro. The logarithm of the density is represented by the color scale shown. The white line shows 
centrifugal barrier for the inflowing gas. The vectors show the poloidal velocity. The plot is for t = 20^0- Simulations 
were performed with resolution Nfi — 150 and Ng — 100. 



through the cross section of the disk. As a result, accretion 
was enhanced inside the main disk as in Model II. Com- 
pared with Model II, the density of counterrotating matter 
is much smaller than that of the "target" disk. For this 
reason the interaction did not lead to dramatic wave for- 
mation. The radial velocity of viscous flow inside the disk 
increased by a factor of ~ 50, and matter flux increased 
by a similar factor. However, the velocity of accretion of 
the main disk remained much smaller than the free-fall 
velocity. 

Note, that during this violent evolution of the outer lay- 
ers of the disk, the main part of the "target" disk was very 
stable, the density did not change appreciably, and the ra- 
dial velocities remained small. 

Figure 18 shows the time evolution of the matter flux 
through the cylindrical surface r — 5i?o- It increases by a 
factor ~ 50 compared with the initial flux. For t/to <J 10 
it increases mainly as a result of surface layer accretion. 
Later, there is also enhanced accretion from the main part 
of the "target" disk. A fraction about ~ 0.2 of the incom- 



ing matter is accreted. 

6. TWO APPLICATIONS 

Here, we briefly consider two applications of our results 
to accretion disks. One is related to galaxies and the sec- 
ond to X-ray pulsars. 

6.1. Galaxies 

In some disk galaxies, counterrotating gas and/or stars 
are observed (e.g., Rubin 1994a, 1994b; and Galletta 
1995). Even though the present work treats Keplerian 
disks, we can nevertheless make estimates based on our hy- 
drodynamic results. The rotation curves of galactic disks 
are flat {v^ =const) over a significant range of distances 
owing in part to a spheroidal distribution {p oc r^^) of 
dark matter. Although this rotation curve is quite differ- 
ent from Keplerian, the analytic viscous accretion flows 
are nearly the same for flat and Keplerian rotation curves 
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Fig. 16. — Model III. Results of simulations of Model III for the case where the inflowing matter rotates in the opposite 
direction to that of the disk. In this case, incoming counterrotating matter interacts strongly with the main disk. The 
resulting low angular momentum gas starts to accrete rapidly along the surfaces of the disk. The times t — 4.5to (top 
panel) and 7to (bottom panel) of evolution are shown. The logarithm of the density is represented by the color scale 
shown. The white lines in both panels show centrifugal barrier for the inflowing gas. The vectors show the poloidal 
velocity. 
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(LC). The reason is that the boundary layer formed be- 
tween the oppositely rotating components, where angular 
momentum is annihilated and fast accretion occurs, does 
not depend strongly on the radial dependence of v^. 

We consider the situation where a counterrotating cloud 
of hydrogen of mass Md = 10^ Mq is captured into an 
equatorial orbit of radius rout = 20kpc. Model III is per- 
tinent to this situation. Interaction between the counter- 
rotating cloud and existing corotating gas (assumed mass 
> Mci) in the galaxy leads to fast boundary layer accre- 
tion along the top and bottom surfaces of the galactic disk 
with average accretion speed ucr ~ 3.25{h/H)Vc, where 
we have adapted the results of §3 to Model III which has 
two boundary layers, and where we have used the circular 
velocity of the galaxy Vc ~ const rather than the Keple- 
rian velocity Vk ■ The time-scale for this fast accretion to 
the center of the galaxy is therefore 

, rout g.. in8 f ro„t \ /200km/s\ 

ic. ~ — ~ 3 X 10 yr [-^) ' (9^) 

where we have assumed h/H = const = 0.1. For the refer- 
ence values of equation (9a) and assuming that the cloud 
mass is accreted in a time ~ lOtcR, the mass accretion 
rate is 

M^^^^0.6MQ/jr, (9b) 

WtCR 

where the factor of two accounts for the fact that half of 
the accreted gas is from the existing disk. Note that the 
accretion time-scale is much less than the Hubble time. 

6.2. X-Ray Pulsars 

Nelson et al. (1997a,b) and Chakrabarty et al. (1997) 
have proposed that the change from spin-up to spin-down 
of X-ray pulsars in binary systems results from a reversal 
of the angular momentum of the wind supplied accreting 
matter. The matter inflowing to the pulsar is assumed 
to be supplied at a constant rate M at a radius rout in 
the equatorial plane of the binary system. Initially, the 
pulsar is assumed to have a conventional, corotating disk. 
Then, due to some disturbance, the sense of rotation of 
the supplied matter reverses so that it is counterrotat- 
ing. Model III is again pertinent to this situation. The 
'new' counterrotating matter encounters the 'old' corotat- 
ing disk. Interaction between the 'new' and 'old' matter 
leads to fast boundary layer accretion with average accre- 
tion speed Ucr ^ 3.25{h/H)VK- The time-scale for this 
fast accretion to reach the pulsar is 

r-j!L.2d(§i)V-!s^)', (10) 

Jo UCR \ M J VlOi^cm/ 

where we have again assumed h/H = const = 0.1. The 
sign of the angular momentum carried by the boundary 
layer flows at the much smaller Alfven radius <C rout, 
which determines the spin-up or spin-down of the pulsar, 
is uncertain from the present work. If this has the sense of 
the counterrotating material supplied at rout, then a rapid 



switch between spin-up and spin-down of the pulsar would 
occur with time scale of order tcR ~ days. 

On a much longer time-scale, inflowing counterrotating 
matter will 'use up' the old corotating disk. This time 
scale is simply the accretion time for a disk rotating in 
one direction. 




A further time interval ^ tss is required for establishment 
of an equilibrium a— disk rotating in the opposite direction 
to the original disk. 

7. CONCLUSIONS 

Time-dependent, axisymmetric hydrodynamic simula- 
tions have been used to study Keplerian accretion disks 
consisting of counterrotating components. The simulation 
code used for the study has a numerical viscosity which 
is calibrated by study of the accretion of a disk rotating 
in one direction. Different grid resolutions were used to 
study the influence of the numerical viscosity. Our simu- 
lation model does not include radiative cooling. Therefore, 
different values of the speciflc heat ratio 7 from 1.01 to 5/3 
were used in order to assess the influence of heating due to 
the numerical viscosity. A small value of 7 — 1 corresponds 
to almost isothermal conditions. 

Different cases were considered. In Model I, the gas 
well above the disk midplane rotates in one direction and 
that well below has the same properties but rotates in the 
opposite direction. In this case there is angular momen- 
tum annihilation in a narrow equatorial boundary layer in 
which matter accretes supcrsonically with a velocity which 
approaches the free-fall velocity. The average accretion 
speed of the disk can be enormously larger than that for 
a conventional a— viscosity disk rotating in one direction. 
For a much lower viscosity (when we took a small part 
of the region and calculated it with a grid resolution of 
200 X 200), the interface between the corotating and coun- 
terrotating components shows significant warping, which 
is probably a type of Kelvin-Helmholtz instability. We 
observed that a large viscosity suppresses this instability. 

In Model II, we considered the case where the inner part 
of the disk corotates while the outer part counterrotates. 
In this case a new equilibrium inner disk forms with a low 
density gap between inner and outer disks. In Model III 
we investigated the case where low-density counterrotat- 
ing matter inflowing from large radial distances encounters 
an existing corotating disk. Friction between the inflow- 
ing matter and the existing disk is found to lead to fast 
boundary layer accretion along the disk surfaces, while in- 
teraction of the disk with counterrotating matter at large 
radii, leads to enhanced accretion in the main body of the 
disk. We observed that the boundary layer accretion is 
a temporary phenomenon, because the interaction of the 
dense disk and low-density counterrotating gas leads to 
heating of this gas. However, the interaction at large radii 
is more steady and leads to continuous enhanced accretion 
in the main disk. 
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These models are pertinent to the formation of counter- 
rotating disks in galaxies, in Active Galactic Nuclei, and in 
X-ray pulsars in binary systems. For galaxies the high ac- 
cretion speed allows coimterrotating gas to be transported 
into the central regions of a galaxy in a time much less 
than the Hubble time. 

It is clear that inclusion of the important role of Kelvin- 
Helmholtz (KH) instabilities in counterrotating disks re- 
quires 3D simulations and this will be the subject of a 
future paper. For example, note that in Model I if the 
flow is treated as a vortex sheet without radial inflow, 
then a local stability analysis k^iJ^ > 1 with k the (r, 0) 
wavevector) indicates unstable KH warping for wavenum- 
bcrs \k^/kr\ < V2{cs/Vk) <C 1 and a maximmn growth 
rate of \kr\cs/2 (LC; Choudhury & Lovelace 1986). To 
account for this non-axisymmetric warping instability of 
the disk, we clearly need 3D simulations. In Model II, 
a local stability analysis of the interface between the in- 
ner and outer disks indicates unstable KH warping for 
\k4,/kz\ < V^Cs/Vk "C 1 and a maximum growth rate 



of \kz\cs/2. Also in this case 3D simulations arc needed. 

A number of further developments are planned. Simula- 
tions explicitly including an a— viscosity are desirable for 
detailed comparison with theory. Explicit inclusion of ra- 
diative cooling is equally important. Also, 3D simulations 
are needed to investigate the Kelvin- Helmholtz instability. 
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Fig. 17. — Model III. Continuation of Figure 16. The top panel is for t — 20to, when the surface layer accretion is strong. 
The bottom panel is for t = 35^0, when the surface accretion is destroyed by the dissipative heating. 
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Fig. 18. — Model III. Matter flux through the cyhndrical surface r = 5i?o in Model III. It has a peak owing to the surface 
accretion, then is almost constant owing to both surface accretion and enhanced accretion inside the main disk. 



